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ABSTRACT 



A new two-dimensional data-adaptive algorithm utilizing the iterative Toeplitz 
approximation method is presented. This algorithm provides a robust and efficient 
means for accurate estimation of 2-D autoregressive parameters representing spatially 
variant data arrays. Its convergence performance is comparable to that of the 2-D 
Recursive Least Squares (RLS) algorithm but is computationally more efficient. The 
savings in computation is realized by reducing the size of the matrix to be inverted 
when solving the AR model normal equations. The ability of the algorithm to accu- 
rately estimate the model parameters using very smcdl data sets and various window- 
ing schemes are evaluated. Spectral estimates are compared for quarter-plane (QP), 
nonsymmetric half-plane (NSHP) and combined-quadrant (CQ) regions of support. 
Additionally, the algorithm is tested in noise cancellation and line enhancement ap- 
plications using image arrays. This algorithm may be implemented for data-adaptive 
image processing or coding and for applications requiring 2-D spectred estimation. 
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I. INTRODUCTION 



A. PARAMETER ESTIMATION OF DATA- ADAPTIVE ARRAYS 

Two-dimensional (2-D) autoregressive (AR) modeling provides a means for the 
spectral estimation of signals received by an array of spatially distributed sensors. 
Additionally, the 2-D AR model paxeuneters may be used to represent the original 
signal in a compact manner resulting in compression of the data to be processed or 
transmitted. This is particularly useful in many image processing applications. The 
traditional leaist squares or Wiener filtering based 2-D AR parameter estimation [Ref. 
1] provide accurate spectral estimates but require the the inversion of a relatively large 
autocorrelation matrix. This is undesirable since matrix inversion is computationally 
expensive. Iterative methods, using a Toeplitz approximation algorithm, have been 
suggested [Ref. 2, 3] as an alternative to the direct Wiener filter solution. The itera- 
tive Toeplitz approximation method has proven to be an effective means to estimate 
the parameters of fixed data arrays without having to invert the correlation matrix. 
This method is particularly useful when data arrays of limited size must be processed. 

The results obtained using iterative methods to obtain the AR parameters of 
fixed data arrays suggest that these methods may be used to some advantage when 
estimating parameters adaptively. Adaptive parameter estimation is necessary in 
many digital signal processing applications where the data is non-stationary. In these 
applications, real-time or near real-time processing is often desirable. This requires 
that any algorithm used must be computationally efficient and be capable of providing 
accurate estimates obtained from a minimum of data. This thesis addresses the 
implementation of the iterative Toeplitz approximation method for 2-D parameter 
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estimation of non-stationary data. It is shown that this algorithm can be efficiently 
implemented to obtain accurate estimates from very small data arrays. 

B. EXPLANATION OF NOTATION 

All rectangular matrices are denoted by a boldface, captiaJ letter, e.g., R. Col- 
umn vectors are designated by a boldface lowercase letter, for example: x. Occasion- 
ally, a vector created temporarily for the development of a mathematical expression 
will be represented by a lowercase, boldface, greek letter. An example of this would 
be the vector a . A special operator matrix will be represented by a calligraphic cap- 
ital letter, such as T. Scalar values are generally represented non-boldface lowercase 
arabic or greek letters, e.g., a or A. Captital letters such as A or P are used to rep- 
resent a variety of values or expressions. These include anything from a polynomial 
to the dimension of a matrix. The symbol T is reserved to indicate the transpose of 
a vector or matrix. 

C. THESIS OUTLINE 

The following describes the organization of the remainder this thesis. Chapter 
II serves as an introduction to the problem of 2-D AR modeling and parameter esti- 
mation. Particular attention is given to the formulation of the normal equations for 
quarter-plane (QP) and nonsymmetric hedf-plane (NSHP) regions of support. This 
chapter concludes with a description of the least squares algorithm and direct in- 
version method of solving the normal equations. This serves as the foundation for 
Chapter III which extends the concept of solving normal equations using the itera- 
tive Toeplitz approximation method to derive parameter estimates referred to as the 
fixed data method. Chapter IV proceeds with the development of the data-adaptive 
algorithm using Toeplitz approximation. Two-dimensional AR spectral estimation 
is discussed in Chapter V, along with a comparison of results using both fixed and 
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adaptive iterative methods for spectral estimation. Additionally, a discussion of the 
use of the data-adaptive algorithm for noise cancellation and line enhancement of 
image arrays is provided in Chapter VI, with experimental results. Conclusions and 
recommendations for future work are found in Chapter VII. 
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II. 2-D AUTOREGRESSIVE MODELING 



A. OVERVIEW 



This procedure assumes a stationary (homogeneous) random process x[ni,n 2 ] 
that is the response of an AR model excited by a white noise input iy[ni,n 2 ] having 
a variance cr^. The AR model as shown in Figure 2.1 may be described as an all pole 



iw(ni,n2) 






R2) 



Figure 2.1: AR model excited by white noise. 



filter with the transfer function 



where A is the region of support over which the parameters a{k\,k 2 ) are non-zero. 
The difference equation for the system that generates x[ni,n 2 ] can be written as 

x[ni,n2] = - -i] + u?[ni,n2]. (2.2) 

(i.i) €>1 

A linear set of equations for the filter coefficients a(i,i) may be formed by multiplying 
both sides of (2.2) by x[ni — l\,n 2 — / 2 ] and computing the the statistical expectation 
of the resulting expression [Ref. 4]. This leads to the following expression called the 
normal equation 



Rr[lu h] = - - i, h - j\ , (2.3) 
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for li,h > 0. The coefficients a{i,j) can be derived from this normal equation. It 
should be noted that the structure of the normal equation depends on the region of 
support A. A region of support may have any shape, but the most commonly used 
are the quarter-plane and nonsymmetric half-plane. 

1. Quarter-Plane Support 

The most straightforward region of support is the quarter-plane. A region 
A is considered to have quarter-plane (QP) support when a{i,j) has non-zero values 
in one quadrant only as shown in Figure 2.2. For this case the normal equation has 




Figure 2.2: The quarter-plane region of support. 



the form 



RAh,h] = e' {ij) ^ (0,0) 

t=0 i=0 



(2.4) 



where /j = 1 , 2, . . . , Px ~ 1 and = 1 , 2, . . . , Pj — 1 with Pi and P 2 being the dimensions 
of A. If it is assumed that a(0, 0) = 1, we may express (2.4) in block matrix forih as 



Ro 


R-i 


R-2 


R 


Ri 


Ro 


R-i 


R 


R 2 


Ri 


Ro 


•• R 



Rpi-l Rpi-2 Rpi-3 



-pl+1 

-Pi +2 
-Pi +3 

Ro 



ao 
ai 
»2 

api-i . 



0 

0 

0 



(2.5) 
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Each block Ra/ of this matrix is given by 



Rm = R^a/ = 



R{M,0) 
R{M, 1) 



R{M,-1) 

R{M,0) 



R{M,-P2 + iy 
R(M,-P2+2) 



R{M,P2-1) R{M,P2-2) 



R(M,0) 



while the coefficients are — [a(M, 0) a(M, 1) a(M, 2) . . . a(M, P 2 ~ 1)]^ and the ^ 
error variance vector is = [<r^ 0 . . . 0]^. The blocks Ra/ along the diagonals of the 
autocorrelation matrix R are equal and the diagonal elements of Ra/ are equal, thus | 
R is Toeplitz block- Toeplitz. 

2. Nonsymmetric Half-Plane Support 

A region A with non-zero a{i,j) as shown in Figure 2.3 is considered to 
have nonsymmetric haif-plane (NSHP) support. The normal equation for this case 




Figure 2.3: The nonsymmetric half-plane region of support. 



must be modified to account for non-zero a(0,j) in the fourth quadrant. This may 
be written as 

LJ+P7-I Fj-lPj-l 

■Rx(/i,/2]= XI o(^J)^T[h,k-j]+ a{i,j)R^[li-i,l2- j] , (2.7) 

i=i t=o j—o 



I 
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where (i,j) ^0, /i = 0, 1, 2, . . . , Pi — 1 and 



. _ J 0, 1,2, . . . , Z /2 + P? ~ for /i — 0 

^ L 2 , L 2 + 1, P 2 + 2, . . . , Zf 2 + P 2 ~ 1, for /i ^ 0 

The dimensions of A are given by Pi and P 2 and L 2 is a negative number as defined 
in Figure 2.3. As in the quarter-plane case, if we let a(0, 0) = 1 we may write (2.7) 
in block matrix form as 



Ro 


R-i 


R-2 


• •• R_p,+r 








-£(0)- 


Ri 


Ro 


R-1 


••• R_P,+2 




ai 




0 


R 2 


Ri 


Ro 


••• R_Pj+3 




®2 


= 


0 


Rpi-i 


RPl-2 


Rpi-3 


• • • Ro 




.ap,_i . 




0 



( 2 . 8 ) 



The blocks Ro, Rm and each have diflFerent forms. The diagonal block 



Ro has the form 



P(0,0) P(0,-1) ••• P(0,-T2-P2 + 1)' 

P(0,1) P(0,0) ••• P(0, -T 2 - P 2 + 2) 

Ro = . . . 

_P(0 , 12 A P 2 — 1) P(0? P 2 d" P 2 — 2) • • • P(0) 0) 

the blocks Rm = R-Im in the first row and column may be written as 



R{M,-L2) P(M,-T2-1) ••• P(M,-L2-P2 + 1)‘ 

P(M,-L2 + 1) P(M,-T2) ••• P(M, -L 2 - P 2 + 2) 

P(M, -T 2 + P 2 - 1) P(M, -Z ,2 + P 2 - 2) ••• P(M,0) 



and the remaining blocks Rm are given by (2.6). 

The model parameter vectors are given by 



ao = [a(0, 0) a(0, 1) a(0, 2) . . . a(0, L 2 + P 2 — 1)]^ 



and 



— [n(-^i -^ 2 ) n(Af, L 2 -f 1) a{M, Z/2 d" 2) . . . a{M, L 2 A P 2 — 1 )]^ 



7 



for M = 1, 2, . . . , Pi — 1. The error variance vector = [cr^ 0 . . . 0]^. Except for the 
upper diagonal block (2.9) and top cind left most blocks (2.10), the NSHP autocorre- 
lation matrix is block-toeplitz with each block being toeplitz as well. Quarter-Plane 
support can be considered a special case of NSHP support with L 2 = 0. 

B. SOLVING THE NORMAL EQUATIONS BY DIRECT INVERSION 
The normal equations also arise in the 2-D linear prediction problem. When 
solving the normal equations, the objective is to find the parameters which mini- 
mize the prediction error [Ref. 5]. It is frequently more convenient to describe the 
formulation of the normal equations in terms of linear prediction. 

In the 2-D least squares problem, the error between the random process x[ni, 712 ] 
and its estimate x[ni,n 2 ] is given by 

e[ni,n2] = x[ni,n2] — x[ni,n2] (2.11) 

or in expanded form as 

e[ni,n2] = a:[ni,n2] + a(t,i)x[ni - t,n2 - j], (i,j) (0,0). (2.12) 

The objective of the least squares method is to minimize the squared errors from a 
particular set of these terms [Ref. 4]. If we let a(0, 0) = 1, we can express (2.12) as 

e[ni,n2] = X^X]°(*’i)®K“*^2-j] (2.13) 

(».j) 

which is compactly represented in vector notation as e = Xa. The error e[ni,n 2 ] 
is computed for each position of the filter, then normal equations are formed by 
multiplying both sides of (2.13) by and applying the orthogonality principle [Ref. 
5). This results in 

Rsl = S (2.14) 
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where R = X^X and the error variance vector S = X^e = 0, . . . , 0]^ with 

£(°) = [cr^, 0, . . . , 0]^. The matrix X may be defined as 





■ 1 




1 1 


1 1 


x = 


xo X' 

1 


= 


Xo Xi • 

1 1 


•• Xp,_i 
1 




1 J 




1 1 


1 J 



and the model parameters may be given by a = [ 1 a']^. The parameter estimates 

are obtained from (2.14) by muliplying both sides of the equation by the inverse of 
R which may be written in expanded form as 

a' = -(X'^X')‘*X'^Xo . (2.16) 

This is referred to as the direct inversion method. 

The preceding discussion assumes that R is Toeplitz block-Toeplitz. This true 
for the case where the autocorrelation method [Ref. 5] is used to form the correlation 
matrix or when the theoretical correlation is known. In practicality, however, pareun- 
eter estimates must be obtained from relatively small sets of finite data. In this case, 
it is often desirable to compute R using the covariance method. This results in a 
block autocorrelation matrix R with non-Toeplitz blocks Rm [Ref. 2]. The rows of 
the matrix X are the elements of the random process covered by the filter mask for 
each point being estimated as the mask is moved over the data. In the covariance 
method the filter mask is not moved past the boundaries of the region of support. 
This means that when QP support is used, X will be of dimension P 1 P 2 x P 1 P 2 where 
Pi and P 2 are the dimensions of the data array accessible to the filter mask. For 
example, using a 3 x 3 (nine element) filter mask to form the normal equations of an 
8x8 data array would result in an X matrix with dimension 36 x 9. It should be 
noted that Pi = P 2 — 6, since two rows and two columns of data cannot be estimated 
by the filter. 

The preceding discussion provides the basis for the subsequent work in this 
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thesis. The next chapter will addresses this same problem using the iterative Toeplitz 
approximation method. 
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III. ITERATIVE TOEPLITZ 
APPROXIMATION ALGORITHM 



The major drawback of the least squares method described in the previous 
chapter comes from the requirement to invert a relatively large autocorrelation ma- 
trix. This is undesireable since matrix inversion is computationally expensive. In 
this chapter we present an iterative method for solving the normal equations which 
does not involve the direct inversion of the correlation matrix to obtain the model 
parameters. This method is then extended to the case where the covariance method 
is used to form the correlation matrix and the normal equations are solved iteratively 
using Toeplitz approximation. 

A. THE ITERATIVE SOLUTION OF THE NORMAL EQUATIONS 

An alternative to direct inversion is to take advantage of the near Toeplitz- 
block- Toeplitz structure of the autocorrelation matrix and successively partition the 
normal equation [Ref. 2]. This partitioning permits the estimation of parameters 
while requiring the inversion of a matrix that is significantly smaller than the original 
autocorrelation matrix. The following discussion develops the iterative solution for 
QP and NSHP regions of support. 

1. Quarter-Plane Support 

In order to iteratively solve the QP normal equations we must begin by 
dividing both sides of (2.5) by <r^. This results in a modified normal equation given 
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by 



Ro 

Ri 

R2 



R-i 

Ro 

Ri 



R-2 

R-i 

Ro 



R-Pi+i 

R-Pi+2 

R-Pi+3 



n 






■5(0)' 




a? 




0 




a^' 


= 


0 




apj_i. 




0 



(3.1) 



.Rpi-1 Rpi-2 Rpi-3 Ro 

where = [1 0 0, . . . , 0]^ and a" = ^a,- for i = 1, 2, . . . , Pi — 1. We now partition 

the normal equation (3.1) as follows 



[Gi 


hi] 






'f 1' 




RoJ 


a" 

L ap-i J 




. 0 . 



(3.2) 



where 



and 





Ro 


R-i 


R-2 


• • • R_Pj +2 ' 




Ri 


Ro 


R-1 


••• R_Pj+3 


Gi = 


R 2 


Ri 


Ro 


R-Pi-^4 




.Rpi-2 


Rpi-3 


Rpi-4 


Ro . 


h[ = 


= [Rpi- 


1 Rpi - 


-2 Rpi- 


3 • • ' Ri ] > 




=\<^ 


_«T 

®i 


a"^ 


a"^ 1^ 

•• apj_2J , 




1 1 = 


[£(0)T 


0^ ... 


. 



(3.3) 



(3.4) 



(3.5) 



(3.6) 



We may now form a set of coupled iterative equations by defining an op- 
erator Pi which is given as 



Pi = 



Gp 0 



0 R^V 

and using this operator to premultiply both sides of (3.2) to yield 



= ®i’h ] -hi a^.,) , 



(3.7) 



(3.8) 



and 



apj_i — Rg^hjai . 



(3.9) 
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Equations (3.8) and (3.9) suggest an iterative solution and can be written as 

= GfM'V 1 -hi , (3.10) 

and 



(3.11) 



a?i. = . 

which provide a means to solve the normal equation (3.1) iteratively. 

It can be noted that the submatrix Gi is nearly the same dimension as 
the original correlation matrix. This means that the inversion of Gi provides little 
advantage, computationally, over the direct inversion method. It is clear, however, 
that further partitioning of Gi will decrease the required complexity. Partitioning of 
Gi yields 



G 2 


h2 


«2 




'* 1 ' 2 ' 


[h^ 


ReJ 


^ P -2 J 




. 0 . 



(3.12) 



where 



and 





Re 


R-i 


R-2 


••• R-Pi+3 




Ri 


Re 


R-1 


••• R_Pj+4 


G 2 = 


R 2 


Ri 


Re 


R-Pi+5 




. Rp, -3 


R'Pj — 4 


Rpi-s 


•• Re - 




= [Rp.- 


2 Rp, . 


-3 Rp,— 


4 • • • Ri ] , 


«2 


= \< 






_//T ]T 
• ^P,-3l 5 




1 2 = 


[£(0)T 


0^ ... 





(3.13) 



(3.14) 



(3.15) 



(3.16) 



As before, both sides of (3.12) are premultiplied by an operator defined by 






•G^-^ 



0 



-1 



(3.17) 
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which results in 

tt2 = 2 ~ ^2 api- 2 ] ^ (3.18) 

and 

= -Ro'h^aa . (3.19) 

Equations (3.18) and (3.19) may be solved iteratively in the same manner as (3.8) 
and (3.9). This solution is expressed as 



= -hjaM, 



(3.20) 



and 

= (3.21) 

This partitioning process may be continued until Gpj_i = R©. This will 
result in a partitioned normal equation given by 



Gpi-i hp,_i] _ [7 

R©Jla^_J-L 0 J’ 



(3.22) 



where hpj_j = Ri, ap,_i = and 7 Pj-i = 

Premultiplication of both sides of (3.22) by an operator will result 

in the iterative solution 



- [R_ia"(*-^^ + . . . + R_p,+,a";i-')] , (3.23) 



and 



a"(*> = -Ro + R-ia''(*-'> + . . . + R_p,+2ay_V^] 

For this case the operator matrix is given by 

T - ° ■ 



(3.24) 



(3.25) 
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At this point in our development, we can combine the preceding iterative 
solutions from the successive levels of partitioning into a set of compact iterative 
summations 

- R,- s' , (3.26) 

1 = 1 

af> = -R^ ■ s’ . (i / i) . (3-27) 

1=0 

where and j = 1, 2, . . . , Pi — 1. The index of 

iteration is k and a" are the P 2 x 1 parameter vectors. The solution of (3.26)-(3.27) 
requires 0{P^P2K) multiplications, where K represents the number of iterations 
to converge to the true parameter values. An additional P 1 P 2 multiplications are 
required to compute the initial pzirameter estimates aj^°^ for j = 1,2, . . . ,P 2 — 1. 
This results in total multplications of O^P^P^K + PiP^)- Depending on the number 
of iterations required to converge to the true parameters, this total compares with the 
algorithms of Wax and Kalaith [Ref. 6] and Akaike [Ref. 7] which require (^(PfPj^) 
operations. When large arrays are used this represents a noticeable savings over direct 
inversion which requires (9(PpPf) multiplications. 

2. Nonsymmetric Half-Plane Support 

The derivation of the iterative NSHP solution follows a procedure of suc- 
cesive partitioning similar to that of QP support. However some differences exist due 
to the asymmetry of the autocorrelation matrix consisting of three different types 
of block matrices. This asymmetry results in a somewhat modified set of recursive 
summations given by 

-W ^ -.(0) _ (3.28) 

t = l 

a"l‘> = -Ra-Ro-'S:'*-'' - R^' s' R,-,a"l‘-'>, (i j), (3.29) 

t=l 
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where and j = 1, 2, . . . , A - 1. The index 

of iteration is k and a" are the P 2 x 1 parameter vectors. As with QP support, The 
solution of (3.28)- (3.29) requires 0{P^PlK) multiplications, where K represents the 
number of iterations to converge to the true parameter values. 

B. THE ITERATIVE SOLUTION USING TOEPLITZ APPROXIMA- 
TION 

The discussion of the previous section assumes an autocorrelation matrix with 
Toeplitz-block-Toeplitz structure. In most cases the autocorrelation matrix must be 
formed from a limited, and possibly very small, amount of sampled data. In this 
situation, it is best to use the covariance method to estimate the correlation matrix 
to reduce the bias in the estimate. However, since this method does not have Toeplitz- 
block-Toeplitz structure, it is necessary to modify the iterative algorithm described 
above. 

1. Forward Iteration 

For first quadrant QP support, the covariance estimate R has the form 



Ro,o 


Ro,i 


Rq,2 


Ro,Pi-i 


Ri,o 


Ri,i 


Ri,2 


Ri,P,_i 


R2,0 


R2,1 


R2,2 


R2,Pi-1 
• • 


.Rpi-1,0 


Rpi-1,1 


RPl-1,2 


••• Rpj_i,Pj_i. 



A proven procedure [Ref. 8] that can be used to form the Toeplitz approx- 
imation is to first average the diagonal blocks Rjj, i = j and then form a Toeplitz 
T matrix from this averaged block Rovg by averaging its diagonal elements. The 
diagonal elements of T can be obtained as follows 

j Pj-t— 1 

~ ~p ^ + iii) (3.31) 

where z = 0, 1, . . . , P 2 — 1- 
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We may now proceed to succesively partition the normal equations in the 
manner outlined in the previous section. The first partitioning leads to the iterative 



equations 



and 



"W _ r'-if 

tti — Gj [7 1 hi apj_j J , 



_ _-R-i 

®Pi-i — “1 

The diagonal block Rp,_i,p,_i may be written as 



Rpi-i,Pi-i = T + Apj_i,p,_] 



(3.32) 



(3.33) 



(3.34) 



where Apj_i,p,_i is the difference between Rp,_i^p,_i ajid the Toeplitz approximation 
T. 

Substituting (3.34) into (3.33) modifies the recursion to give 

a S‘* = G.-’ h , - h. a"S‘.-’>] , (3.35) 

and 

apfii = -T-'Ap,_i.p,_ia;f_-'^ - . (3.36) 

The final expressions resulting from the successive partitioning are then 

a p, -1 = G?,'-. h n - hp, a;'"-’') , (3.37) 

and 

a"<''> = -Ri,ihp,_ial^-V. (3.38) 

where apj_i = ao,Gpj_i = Ro,o» and hp,_i = Ro,i- The diagonals Hjj may be 
rewritten as 

Rjj = T + Ajj , (3.39) 
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Now, substituting (3.39) into (3.38) results in 

— hpj_i , (3.40) 

and 

^ - T-^hfapJiV^ • (3-41) 

Finally, combining the succesive iterative equations results in our solution to the 
normal equations using the forward iterative algorithm for the QP case, which tahes 
the form 

= ai'l”' - T-'A„, „<“-■> - T‘' (3.42) 

t = l 

a''<‘) = -T-’ - T"' (3.43) 

•=0 

where = T“^Rj,oa^^^°^ and j = 1,2, . . . , Pi — 1. The forward 

iteration will provide parameters for the first quadrant spectral estimate. 

2. Backward Iteration 

A filter mask may have quarter-plane support in any two quadrants. A sec- 
ond quadrant estimate exists for the quadrant adjacent to either side of the quadrant 
designated to be the first quadrant. The Hermitian symmetry and Toeplitz-block- 
Toeplitz nature of the autocorrelation matrix causes the estimates produced from 
diagonally adjacent quadrants to be identical [Ref. 9]. 

The second quadrant AR parameters b{i,j) may be estimated from the nor- 
mal equations using the backward iterative Toeplitz approximation. When 6(0, 0) = 1 
the second quadrant normal equation is given by 



1 

b 

K> 

1 




bp,_i ■ 




■ 0 


B-1,0 R-1,1 R-1,2 R-i,Pi-i 




bp, -2 




0 


R-2,0 R2,1 R2,2 *•* R-2,Pj-1 




bp, -3 


= 


0 


-Rpi-1,0 Rpi-1,1 Rpi-1,2 ••• RPi-l,Pi-l- 




. bo . 




£(0) 
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where R,j = Rj,-, = [100,..., 0]^ and 



bl = (l/^")[fc(^, 0), b{k, 1), b{k, 2), . . . , b{k. Pi - 1)] . 



It may be noted that the second quadrant autocorrelation matrix is identical to that 
of the first quadrant. However, the parameter estimates b{i,j) are not generally the 
same. The estimation of the backward parameters is computed iteratively in the same 
manner the forward parameters, but, the backward method differs in the way that 
the normal equation is partitioned. In this case the partitioning begins at the upper 
left diagonal block and continues until the G matrix consists of only the lower right 
diagonal block. The Toeplitz approximation T is identical to that of first quadrant 
support. 

The combined backward iterative equations, which are similar to those for 
the forward iteration, can be written as 



t = l 


(3.45) 


Pl-1 


53 Rpi-i-j.Pi-i-.b>j_V_, 

t=0 


(3.46) 



where bf = T-iRp^.i^bf and j = 1, 2, . . . , Pi - 1. The first and 

second quadrant parameters are used to find the combined- quadrant spectral estimate, 
which is discussed in Chapter V. 

3. Iterative Method for NSHP Support 

The iterative Toeplitz approximation method may be used to estimate the 
parameters of arrays with nonsymmetric half-plane support. As shown previously the 
NSHP autocorrelation matrix ha^ an asymmetric structure which can be expressed 
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(3.47) 



^,0 


Ro,i 


Ro,2 


• * • He, Pi -1 


Ri,o 


Ri,i 


Ri,2 


Ri,Pi-i 


6-2,0 


Rz.i 


R2,2 


R2,Pj_i 


-Rpi-1,0 


Rpi-1,1 


RPl-1,2 


••• Rpi-i,Pi-i. 



The dimensions of the top most and left most blocks are reduced by L 2 rows and 
columns, respectively. As in the previous cases the iterative process begins with 
forming the Toeplitz approximation of Rawp. In the NSHP case, however, the upper 
left main diagonal block is not included in the average becuase it does not have 
the same dimension as the other mcdn diagonal blocks. Thus, it is necessary to use 
the elements of the upper left main diagonal block to create a separate Toeplitz 
approximation designated as T. 

The normal equation is successively partitioned as described in the previous 
sections. At each stage of the partitioning, the blocks along the main diagonal are 
expressed as the sum of the Toeplitz approximation and a difference matrix. The 
iterative summations are given by 



- T-> y: 
1=1 


(3.48) 


Pi-i 


-1) _ T-i Y: 
«=0 


(3.49) 



where = T and j = 1,2, . . . , Pi — 1. 



C. SUMMARY 

It is clear that iterative methods have a computational advantage over the di- 
rect inversion method. This has been shown for the Toeplitz-block- Toeplitz case and 
extended to the case where Toeplitz approximation is used to compensate for the 
non-Toeplitz nature of autocorrelation matrices formed using the covariance method. 
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Generally, the matrix that must be inverted using the iterative method has the di- 
mension of a single block in the block matrix. This corresponds to the filter order or 
mask size. The mask size will vary with the order necessary for accurate AR mod- 
eling of the sampled data, but filter masks that range from 3 x 3 to 6 x 6 will be 
sufficient for most applications. Although other methods may exhibit faster conver- 
gence to the true parameters, they do so with greater computational complexity. The 
iterative methods provide an approximation much sooner. Additionally, the iterative 
methods are guaranteed to converge when the block matrix is symmetric and positive 
definite [Ref. 2]. This property is important when extending the iterative Toeplitz 
approximation algorithm to the data-adaptive case which is introduced in the next 
chapter. 
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IV. DATA-ADAPTIVE ITERATIVE TOEPLITZ 
APPROXIMATION ALGORITHM 



A. OVERVIEW 

The necessity for adaptive filtering techniques arises when it is desired to process 
signals that result from an environment of unknown statistics [Ref. 10]. There exists a 
broad class of problems that fall into this category, these include such diverse fields as 
sonar, radar, image processing, seismology and communications. In general, adaptive 
filters provide a significant improvement in performance over fixed coefficient filters 
designed to operate on data with unknown statistics. Additionally, the use of adaptive 
filters makes possible new signal processing capabilties that would not be available 
otherwise. 

One distinct advantage associated with adaptive filters pertains to their ability 
process data on-line. This is desirable for many applications such as autoregressive 
spectrum analysis, detection of a signal in noise, adaptive noise cancellation and line 
enhancement. 

Two-dimensional adaptive filters are used to process data obtained from spa- 
tially variant data arrays. The most successful algorithms currently implemented for 
this purpose include the 2-D least mean squaxe (LMS) and recursive least squares 
(RLS) algorithms. Both of these algorithms have been derived from the 2-D Wiener 
filter and method of least squares [Ref. 11). With this in mind, it follows that the 
iterative method for fixed data may be successfully extended to operate on spatially 
variant arrays. The remainder of this chapter describes the development of the data- 
adaptive iterative Toeplitz approximation algorithms and some considerations that 
apply when using this algorithm. 
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B. THE DATA-ADAPTIVE ALGORITHM 



The data-adaptive Toeplitz approximation algorithm is based on the same suc- 
cesive partitioning described for the fixed data case. The algorithm becomes adaptive 
in nature when the autocorrelation matrix Rn is updated for each shift n of the filter 
mask. This yields a continuously updated set of AR parameter estimates. Compu- 
tation of AR parameter estimates begins with the first iteration [Ref. 12]. As in the 
fixed data case, only the approximated Toeplitz matrix T is inverted. Computation 
time is shortened further by the fact that R„ is computed only from data covered 
by the filter mask, which amounts to taking the outer product of a P 1 P 2 x 1 vector 
where P 1 P 2 is equal to the number of filter coefficients. The data-adaptive form of 
the combined iterative equations can be written as 

«;<“> = a;"”' - T;'A„, - T;> (4.1) 

t=l 

a"(”) = -T;' - T;> s' (4.2) 

i=0 

ijtj 

where j = 1,2, . . . , Pi — 1. These can be 

compared to (3.42) and (3.43) for the fixed data case. The initial parameter estimates 
aio" are determined for the initial point being estimated and all subsequent parameters 
are a function of the parameters computed at the previous meisk positions. The 
backward equations, used for the second quadrant estimates, can be written as 




b(0) _ 




t = l 




Pl-1 



-T;'A„.p,_._,,p,...ib;"-"-T;' S 



•=0 



(4.3) 

(4.4) 



where b<“> = b„f = T„-‘R<,p,.,,^b'‘'> and j = 1,2,. . . ,P, - 1. 

It should be noted that the only difference between the data-adaptive equations 
and the fixed data equations is the the index of iteration which is n instead of k 
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and the manner in which the autocorrelation matrix is formulated. In the data 
adaptive case, the iteration proceeds with each shift of the filter mask or update of 
the autocorrelation matrix which has the recursive form 

= R„_i + x„x^ . (4.5) 

The P 1 P 2 X 1 vector x„ is obtained by arranging the 2-D data covered by the filter 
mask at each shift n = (ni,n 2 ) in a vector form. The diagonal blocks of Rn are used 
to form T„ in the same manner as the fixed data case. 

As with the RLS algorithm and other recursive implementations of the method 
of least squares, it is necessary to introduce a weighting or forgetting factor when 
computing Rn [Ref. 10]. The forgetting factor is a scalar value that may be designated 
as 0{n) where n is the iteration number of the point being estimated. The weighting 
factor P{n) has the property that 

0 < y3(n) < 1, n = l,2, ...,p, (4.6) 

where p is the total number of points being estimated. The purpose of ^{n) is to 
ensure that data in the distant past is ‘forgotten’ which will make it possible for 
the filter to adapt to statistical variations of the observed data when operating in a 
non-homogeneous environment. A commonly used form of the weighting factor is the 
exponential weighting factor which is defined as 

/?(n) = AP~”, n = l,2,...,p, (4.7) 

where A is positive constant with a value less than 1. The quantity 1/(1 — A) can be 
considered as an approximate measure of the memory of the algorithm. When A = 1 
we have the ordinary method of least squares which corresponds to infinite memory. 
By this reasoning, we may introduce a method of exponentially weighted least squares, 
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where ^{n) is used as a weighting factor in the expression of the performance measure 
[Ref. 10] 

{(n) = X;A»-”|e(n)p, (4.8) 

n=l 

where e(n) is the error defined by (2.12) at mask position (ni,n 2 )- 

To see how the forgetting factor is introduced in the recursion we must refer to 
the original definition of the normal equation 

Ra = S , (4.9) 



where a is the vector of optimum parameters which results in the minimum value of 
the performance measure ^(n). The P 1 P 2 x P 1 P 2 autocorrelation matrix in this case 
is defined by 

R, = x: AP-"x„x^ (4.10) 

n=l 

We may now isolate the term of (4.10) that corresponds to n = p from the rest of the 
summation on the right side of the expression to obtain 

Rn = A 

The term inside the brackets on the right side of (4.11) is actually the sum of the pre- 
viously computed autocorrelation matrices R(n — 1), therefore we have the recursive 
relationship 

R„ = ARn-i + x„x^ . (4.12) 

The affect of /3(n) is apparent when (4.12) is expanded to give 



X^A^ ^ "Xn_lX^_l 
n=l 



+ . 



(4.11) 



Ri = xixf 

R2 = Axixf -f X2X2 

R3 = A^xixf + AX2X2 + X3X3 



R„ = A””^xixf -I- A"“^X 2 X^-f ... + x„x^ (4.13) 
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which clearly shows the decreasing contribution made by past data vectors when 
A < 1. 

The choice of a value for A is dependent on the statistical nature of the data 
being processed. Generally, if the data is random or uncorrelated it is undesirable to 
use the past data to compute the current estimate, therefore a small value for A would 
be the appropriate choice. Conversely, a value of A close to 1 would be desirable for 
highly correlated data. 

Since the statistical nature of the data is often not known a priori when using 
data-adaptive filters, it is difficult to choose an optimum value for the forgetting 
factor. A further modification of the recursive equation (4.12) has been implemented 
for the data-adaptive iterative algorithm to address this problem. This modification 
involves the application of a weighting factor to the current update of Rn- This 
weighting factor depends on the index of iteration n such that its value approaches 
1 as n — > oo. This biasing scheme causes the influence of the present data on the 
parameter estimates to increase with the number of iterations. The modified recursion 
is given by 

Rn = ARn-l + (~“) • (4-14) 

The weighting factor on the second term in (4.14) is added to provide a more gradual 
shift of influence to recent updates than would be possible with only the A term 
present. This factor is not as important when a highly correlated signal is processed. 

C. SUMMARY 

The adaptive Toeplitz approximation algorithm can be used to estimate param- 
eters for the same regions of support as in the fixed data case. Additionally, various 
schemes for pre-windowing or post-windowing the data may be employed to improve 
the estimates. Other factors affecting the performance of the algorithm include the 
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directional scheme for moving the filter over the data array and the criteria used to 
determine convergence of the parameter estimates. This algorithm heis been used for 
spectral estimation, image noise cancellation, and line enhancement. Experimental 
results involving these applications are presented in the next two chapters. 
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V. 2-D AR SPECTRAL ESTIMATION 



A. OVERVIEW 

One application where the use of two-dimensional autoregressive modeling has 
proven to be paxticularly advantageous is that of spectral estimation. The 2-D AR 
model Wcis derived in Chapter II. This model is used to find the 2-D AR power 
spectral estimate Px{u>i,(jJ 2 ) which is given by [Ref. 4] 

Px(u;,,a;2) = |if(u;„a;2)|2p^ , (5.1) 

where is the power spectral density of the input and is the transfer 

function of the 2-D AR model. The input is white noise with a constant power 
spectrum of amplitude Therefore, we can write (5.1) as 

2 

+ ■ (5.2) 

The key to finding the 2-D AR power spectral estimate is determining the the pa- 
rameters a{k\,k 2 )- The next section provides results of experimentation with the 
data-adaptive iterative algorithm as used to find the parameter estimates of a 2-D 
signal consisting of a sum of sinusoids in additive white noise. Results using the di- 
rect inversion method and fixed data iterative Toeplitz approximation are provided 
for comparison. 

B. EXPERIMENTAL RESULTS 

The performance of the algorithms discussed in the previous chapters is com- 
pared by examining the estimated spectral densities computed from the parameters 
produced by each method. Each algorithm has been used to estimate the parameters 
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of the sinusoidal input array generated by 

x(ni,n2) = cos(27r/iini + 27r/i2n2) 

+A 2 cos(27r/2ini + 27: f 22 ^ 2 ) + 'w{ni , n2) (5.3) 

where amplitudes Ai = A 2 = y/2 and fij axe normalized frequencies in the range 
(0 ^ fij ^ 0.5). The frequencies chosen are /n = /12 = 0.167 and /21 = f 22 = 0.333. 
Additionally, the data-adaptive algorithm is tested with off-set frequencies to evaluate 
its peformance for that case. The noise term w{ni,n 2 ) is zero mean and gaussian with 
the variance cr^ scaled give a desired signal-to-noise ratio (SNR). The SNR (in dB) 
is defined as [Ref. 4] 

SNR=101og,„l;^. (S-l) 

where N is the number of sinusoids present and A is the amplitude term. The SNR 
is varied by holding At constant and adjusting the value of cr^. The algorithm’s 
performance using a 3 x 3 filter mask is evaluated for various sizes of input arrays 
and the common regions of support, including combined-quadrant (CQ). The results 
obtained using the data-adaptive algorithm are provided for the pre-windowed and 
covariance methods. The data-adaptive algorithm is tested with SNR’s of 10 dB and 
0 dB. It should be noted that the surface plots in all figures have been rotated 90 
degrees to show the separation of the spectral peaks more clearly. Contour plots are 
provided to better judge the accuracy of the estimates. The axes of the contour plots 
range from 0-30 which is taken from a 32 x 32 frequency domain array representing a 
quadrant of a 64 x 64 point 2-D FFT output. The 32 x 32 array covers the frequency 
range from 0 to tt. 

1. Estimates by Direct Inversion 

The first case examined is that involving direct inversion of the autocorre- 
lation matrix to find the least squares AR parameters from the normal equations. An 
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8x8 input array was used with 1" quadrant quarter-plane support. Tlie resulting 
spectral density estimate is shown in Figure 5.1. This plot shows spectral peaks at 





Figure 5.1: Spectral estimate using direct matrix inversion; SNR=10 dB. 



21.3 and 10.7 on the FI and F2 axis respectively, which correspond to the normalized 
frquencies fn = /ij = 0.167 and /n = /jj = 0.333 of the input data. The true 
locations of the frequencies are indicated by crosses. 

2. Spectral Estimates using the Iterative Method for Fixed Data 
a. Quarter-Plane Support 

Figure 5.2 shows the spectral density of x-(ni,n 2 ) cis estimated using 
the fixed data iterative method with first quadrant QP support. As in the previous 
case the frequencies estimated are very close to the true frequencies. It is interesting 
to note however, that contrary to what might be expected, the quality of the estimate 
appears to degrade as the number of iterations is increased. The estimate produced 
after one iteration is clearly superior to the estimate obtained after ten iterations. 
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This plieiionieiion is consistent witli tlicit experienced in the 1— D case wlieie the 
iterative solution would begin to diverge after a certain number of iterations in cases 
with small data sets [Ref. 12]. 
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Figure 5.2: Fixed Data, 1st Quadrant QP Support, SNR=10 dB. 

Results using second quadrant QP support proved to be generally 
better for this placement of the sinusoids than those obtained from the first quadrant. 
An example is provided in Figure 5.3. 

b. Combined Quadrant Support 

Spectral density estimates using QP support demonstrate a tendency 
to distort elliptically. This distortion is evident in the first quadrant QP case shown 
in Figure 5.2. A method to compensate for this elongation of the spectral peaks 
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Figure 5.3: Fixed Data, 2nd Quadrant QP Support, SNR=10 dB. 



32 



involves combining the first and second quadrant estimates [Ref. 13) to give 

2^2 



P. = 



i’,' + Ai'l 



(5.5) 



The plots in Figure 5.4 indicate that this method yields accurate estimates. However, 
despite some improvement, the elliptical elongation at the base of the spectral peaks 
resulting from the influence of the first quadrant estimate is still quite noticeable. 




10 licrailona 1 Iteration 





Figure 5.4: Fixed Data, Combined Quadrant Support, SNR=10 dB. 



3. Spectral Estimates using the Data-Adaptive Iterative Method 
a. Pre-Windowed Input Arrays 

The input arrays were pre-windowed by adding Pj — 1 rows and 
columns of zeros to the top and left side of the array, where P\ is the dimension 
of the filter mask. The purpose of pre-windowing the data is to compute the pa- 
rameter estimates using all the input data. Figure 5.5 shows the spectral estimates 
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for 1*', 2"^^ and CQ support with an 8 x 8 input array. In Figure 5.6 we present the 




Figure 5.5:' 8 x 8 input array (pre-windowed), SNR=10 dB. 



spectral estimates when only 16 data points (4x4 array) were available. The spectral 
estimates obtained from the 8x8 input array were quite acceptable, but some bias 
is evident in the first quadrcuit estimate which is carried over to the CQ estimate. 
The estimates obtained from the 4x4 array demonstrate the value of computing the 
CQ estimate. In this case the true signal frequencies were unrecognizable until the 
combined quadrant estimate was obtained. 

A final test of the pre-windowed algorithm consisted of running the 
filter mask only half way through an 8 x 8 input. The estimates are computed from 
the upper 4x8 section of the observed data. The result of this experiment is given in 
Figure 5.7. These results are similar to that of the 4x4 Ccise and serve to reinforce the 
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Figure 5.6: 4x4 input array (pre- windowed), SNR=10 dB. 
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Figure 5.7: 4x8 portion of the input array (pre-windowed), SNR=10 dB. 
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importance of the Combined Quadrant technique for 2-D AR spectrum estimation, 
b. Estimates using the Covariance Method 

As explained in section B of Chapter II the covariance method uses 
only the observed data. As a result, Pi — 1 rows and columns of the input array are 
not available when computing the parameter estimates. However, the experimental 
results obtained in this work indicate that the covariance method achieves better 
spectral peak resolution. Additionally, less bicis is observed in the spectral estimates. 
This characteristic of the covariance method is evident in Fig 5.8. In this Ccise the 




Figure 5.8: 8x8 input array (Covariance Method), SNR=10 dB. 

best spectral estimate is clearly provided by the Combined Quadrant method. The 
bias present in the pre-windowed first quadrant support example haa been completely 
eliminated by using the covariance method. 
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The next example demonstrates a unique characteristic of the data- 
adaptive iterative Toeplitz approximation algorithm. In this example the spectral 
estimate is computed using a 3 x 3 filter mask with a 4 x 4 input array. This means 
that only four elements of the input data are available to the filter to determine 
the parameter estimates. This result is of particular interest because the correlation 
matrix is actually rank deficient, and the problem cannot be solved by direct inversion. 
Nevertheless, the iterative Toeplitz approximation method does provide a means to 
compute the estimate. The results are given in Figure 5.9. Although the spectral 
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Figure 5.9: 4x4 input array (Covariance Method), SNR=10 dB. 

estimates for the first quadrant and Combined Quadrants axe somewhat weak, the 
second quadrant estimate exhibits remarkable resolution and accuracy. 

As a further test of the algorithm, a signal with offset frequencies 
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fn = 0.250, /i 2 = 0.125, / 2 i = 0.300, and /22 = 0.400 is generated and processed. 
The spectral estimates resulting from the computed model parameters are plotted 
in Figure 5.10. For this case, the best spectral estimate was obtained using first 




Figure 5.10: 8x8 input array, offset frequencies, (Covariance Method), 
SNR=10 dB. 

quadrant suppport. 

The spectral estimates for a signal-to-noise ratio of 0 dB is ploted 
in Figure 5.11. The accuracy of the estimates is noticeably degraded in this case. 
However, a reasonable indication of the signaj frequency content can be seen when 
the combined-quadrant method is used. 

The final c<ise using NSHP support is included for completeness. The 
filter mask had dimensions Fj = F 2 = 4 and = —2. The parameters were estimated 
for a. 16 X 16 input array. The NSHP spectral estimate is plotted in Figure 5.12. The 
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Figure 5.11: 8x8 input array (Covariance Method), SNR=OdB. 




Figure 5.12: 16 x 16 input array (Covariance Method), SNR=10 dB; non- 
synimetric half-plane support. 
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resulting spectral estimate shown in Figure 5.12 is quite accurate, but this comes 
at the cost of greater complexity for the filter. NSHP support may be used as an 
alternative to combined- quadrant support in some cases. 

C. SUMMARY 

The data-adaptive iterative Toeplitz approximation algorithm has proven to 
be a viable means for 2-D AR spectrum estimation. Numerous factors must be 
considered when using this algorithm. In particular, the filter mask size should be 
chosen to correspond to the model order of the signal being processed. Although 
our overall experience has found that that CQ support produced the best results, 
no region of support seemed to have a clear advantage over the others for all cases. 
This indicates that the optimum region of support must be evaluated for the specific 
application that the algorithm is implemented for. 
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VI. RESTORATION OF IMAGES 



A. OVERVIEW 

The characteristics of typical images vary greatly from one region of the image 
to the next. For example, a region of an image containing a crowd or a building may 
have detailed variations in intensity, while another region representing the sky in the 
background will be essentially uniform in intensity. Two-dimensional data-adaptive 
filtering is particularly suited for processing data of this nature. The idea of axiapting 
the processing to the local characteristics of an image is advantageous for many image 
processing applications, including image enhancement and restoration. 

There are two approaches to adaptive signal processing. The first approach 
involves adapting the iilter output for each pixel. This is known as pixel-by-pixel 
processing [Ref. 4]. In this scheme, the processing method is based on the local 
characteristics of the image, degradation, and any other pertinent information con- 
tained in the pixel’s neighborhood region. This approach offers the greatest degree 
of flexiblity to the adaptive process, but has the highest computational complexity. 

When a more computationally efficient method is desired, a second approach 
known as subimage-hy-subimage or block-by-block processing can be used [Ref. 4]. In 
this approach, the image array is divided into subimages and space-invariant tech- 
niques are used to process the data. This method can result in some discontinuities 
being present in the processed image. The size of the subimages directly affects the 
quality of the output image. The block-by-block method provides a faster means of 
processing, however, and may be desirable in some applications. 

In the next two sections of this chapter the data-adaptive Toeplitz approxima- 
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tion algorithm is implemented in image noise cancellation and line enhancer modes. 
For both of these applications the algorithm has been modified to perform pixel-by- 
pixel processing. 

B. IMAGE NOISE CANCELLATION 

The usual method of finding the estimate of a signal in noise is to pass the 
corrupted signal through a system that serves to suppress the noise while leaving 
the desired signal relatively unchanged. The noise canceler developed by Widrow 
[Ref. 14] is an example of such a system. Adaptive noise cancellation is a variation 
of optimal filtering that can be used for image restoration to some advantage. In 
particular, when a recieved image is obscured by noise or another image transmitted 
from another location as shown in Figure 6.1, an adaptive noise canceler can be 
employed to extract the desired image from the composite signal. The noise canceler 



unwanted transmission 




Figure 6.1: A possible scenario for application of a noise canceler as used 
for image restoration. 

makes use of an auxiliary or reference input derived from one or more sensors located 
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at points in the noise field where the desired signal is weak or undetectable. This 
input is filtered and subtracted from the primary input containing both desired and 
unwanted signals. A block diagram of the noise canceler is shown in Figure 6.2. The 
reference signal and the undesired part of the primary input are correlated. For this 
reason, the unwanted signal or noise is eliminated or attenuated by cancellation. 




Figure 6.2: Adaptive Noise Canceler Block Diagram. 



To evaluate the performance of the data-adaptive iterative Toeplitz approxi- 
mation algorithm used in a noise canceler configuration, a 256 x 256 composite im- 
age array weis created from two distinct images. This corrupted image is shown in 
Figure 6.3. The algorithm described in Chapter IV wais modified to include a cross- 
correlation term given by 

r„ = Ar„_i -f «^nX„ (6.1) 

where d„ is the element of the composite signal corresponding to the element of 
the reference signal being processed by the filter. This cross-correlation term is used 
to form the initial parameter estimate for each pixel being processed. 
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Figure 6.3: Corrupted image array provided as primary input to noise 
canceler. 

The resulting error e(ni,n 2 ) computed at the output of the system becomes the 
restored image. This relationship is given by 

e(ni,n2) = [d(ni,n2) + it(ni,n2)] - u(ni,n2) , (6.2) 

where d is the desired image, u is the unwanted image or noise and u is the filtered 
estimate of u. 

In the simulation the noise signal added to the desired image is scaled and 
lowpass filtered to simulate the effect of having the primary and reference sensors 
physically separated. The output of the noise canceler is shown in Figure 6.4. The 
unwanted image has been noticeably attenuated and the desired image is clearly 
visible. The original images are provided in Figure 6.5 as a reference for evaluation 
of the algorithm’s performance. 

C. ADAPTIVE LINE ENHANCER 

A special case of adaptive noise canceling occurs when only the signal that is 
corrupted by noise is available. In this case the recieved signal is delayed and used as 
the reference signal. The reference signal may be expressed as t/(ni, ri 2 ) = ar(ni — 6 ,n 2 ), 
where 6 is the amount of spatial delay used. The block diagram for this system is 
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Figure 6.4: Processed image array after noise cancellation. 




Figure 6.5: The original images - undesired and desired. 
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shown in Figure 6.6. The main function of the delay parameter 8 is to remove any 




Figure 6.6: The adaptive line enhancer block diagram. 

correlation that may exist between the noise component in the original input signal 
and the noise component of the delayed adaptive filter input [Ref. 10]. The filter will 
respond by cancelling any components of the main signal x(ni,n 2 ) that are in any 
way correlated with the secondary signal y(ni,n 2 ) [Ref. 11]. The remaining noise 
component at the output of the filter is then subtracted from main signal resulting 
in the removal of the additive noise in the signal. In general, the delay 8 should 
be chosen such that it is approximately half the length of the correlation sequence 
[Ref. 11], The input image and output of the adaptive line enhancer are shown in 
Figure 6.7. 

D. SUMMARY 

In this chapter, it has been demonstrated that the data adaptive Toeplitz ap- 
proximation algorithm may be used in image processing applications. While pixel-by- 
pixel jiiocessing was used in the above examples, it may be desirable to combine this 
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Figure 6.7: Input and output of the data-adaptive line enhancer. 



approach with block-by-block processing to obtain better results. This would entail 
dividing the image into smaller blocks before processing and then recombining the 
processed blocks. Additionally, reprocessing of the array estimates may be a viable 
means of improving image quality. 
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VII. CONCLUSIONS 



The data-adaptive iterative Toeplitz approximation algorithm is readily im- 
plemented from the iterative methods previously used for fixed data-arrays. The 
experimental results indicate that this algorithm is well-suited for a wide variety of 
2-D signal processing applications. Practically speaking, the iterative method can 
be successfully implemented for any application where the use of Wiener filter based 
adaptive filters has been successful. While more extensive evaluation is necessary 
before this method can be said to have a clear advantage over other methods, there 
is enough evidence of its capabilities to warrant further investigation in this area. 

A. PERFORMANCE EVALUATION SUMMARY 

As stated earlier, autoregressive model parameters may used to estimate the 
spectral content of 2-D arrays with high resolution. In experiments the performance 
of the data-adaptive iterative algorithm matched that of the direct inversion method, 
and in some cases, improved upon the performance of the fixed data iterative method. 
This improvement may be due to the way in which the iteration is carried out in the 
adaptive case versus the fixed case. In the fixed case the iteration is carried out using 
a correlation matrix formed from all the data at once. In the data-adaptive case the 
correlation matrix is formed recursively with new data incorporated at each itera- 
tion. Of particular note, was the ability of the data-adaptive algorithm to estimate 
frequencies using very small data sets and its ability to produce an estimate even 
when the correlation matrix is not of full rank . Additionally, its performance in a 
high noise environment was noteworthy. 

The algorithm also proved to be a viable method for data-adaptive image 
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restoration. The nature of this application differs greatly from that of spectral estima- 
tion in that the data sets are much larger and generally less correlated. Encouraging 
results were obtained for the data-adaptive noise canceler and line enhancer with only 
minor modifications of the same algorithm that was used for the spectrum estimates. 

B. RECOMMENDATIONS FOR FUTURE STUDY 

The primary objective of this thesis was to develop the data-adaptive iterative 
Toeplitz algorithm and obtain experimental results to evaluate its potential for use in 
applications that involve 2-D AR parameter estimation. There are many aspects re- 
garding this work that require more in-depth study and several applications that may 
benefit in the use of this algorithm. In particular, more detailed analysis of model 
orders and their relation to the algorithm’s ability to form spectral estimates from 
very small input arrays needs to be carried out. The use of singular value decompo- 
sition (SVD) and eigenvalue analysis of the observed data may provide some insights 
regarding the performance of the iterative algorithm. Also, further investigation as to 
why one region of support (ROS) provides greater resolution than another, may yield 
information that would permit the choice of an optimum ROS prior to processing 
the data. Further analysis for choosing values of the forgetting factor A and how it 
relates to the statistical nature of the data, is desirable as well. This may lead to a 
scheme for on-line adjustment of A, which would be particularly advantageous when 
processing images. Other future work could explore the use of different data win- 
dowing schemes and data-ordering schemes for passing the filter mask over the data. 
One possible method that has not been investigated is an expanding square starting 
at one corner of the array. This method seems well suited for block-by-block image 
processing. Finally, the data-adaptive iterative Toeplitz approximation algorithm has 
the potential for use in non-linear applications, such as the identification of Volterra 
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